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The weakly nonlinear regime of a viscoelastic Navier-Stokes fluid is investigated. For 
the purely hydrodynamic case, it is known that large-scale perturbations tend to the 
minima of a Ginzburg-Landau free-energy functional with a double-well (fourth-order) 
potential. The dynamics of the relaxation process is ruled by a one-dimensional Cahn- 
Hilliard equation that dictates the hyperbolic tangent profiles of kink-antikink structures 
and their mutual interactions. For the viscoelastic case, we found that the dynamics still 
admits a formulation in terms of a Ginzburg-Landau free-energy functional. For suffi- 
ciently small elasticities, the phenomenology is very similar to the purely hydrodynamic 
case: the free-energy functional is still a fourth-order potential and slightly perturbed 
kink-antikink structures hold. For sufficiently large elasticities, a critical point sets in: 
the fourth-order term changes sign and the next-order nonlinearity must be taken into 
account. Despite the double- well structure of the potential, the one-dimensional nature 
of the problem makes the dynamics sensitive to the details of the potential. We analysed 
the interactions among these generalized kink-antikink structures, demonstrating their 
role in a new, elastic instability. Finally, consequences for the problem of polymer drag 
reduction are presented. 



1. Introduction 

The derivation of coarse-grained equations of motion, averaging out microscopic de- 
grees of freedom and retaining only those features relevant to the process of interest, is a 
major goal in many different scientific domains. A first classical example is the dynamics 
of celestial bodies, the physical problem which motivated the introduction of asymptotic 
techniques to systematically average over rapidly rotating, angular degrees of freedom. 
More recently, many interesting phenomena in biological contexts (e.g. related to do- 
main formation in lipid membrane, bilayer fusion, and cooperative motions associated 
with phase changes) have been found to occur on times and length scales much larger 
than the typical times and scales where t he classical molecular-dynamics methods are ap- 
plicable l|Vat,t,ulainen & KarttuneiJ l2005). To reach those larger length-scales, one resorts 
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to co arse-grained models that employ effective interaction potentials ijKarttunen et all 
l2004 . 

Another relevant example of coarse-grained model comes from climatology. The current 
numerical models for prediction of weather and climate involve general circulation mod- 
els. They consist of coupled, nonlinear partial differential equations, discretized in space 
and time for the purpose of numerical simulations. The current generation of supercom- 
puters supports mesh spacing of the order of 200 km for short-term climate simulations. 
However, many important physical processes occur on smaller scales (e.g. the cloud cover 
in the boundary layer) and they significantly affect the large-scale dynamics of resolved 
fields. A powerful way to incorp orate the unresolved d ynamics is provided by suitable 
coarse-grained stochastic models l|Khouider et QiJf2OO.il) . 

Finally, in the framework of phase-ordering kinetics, the concept of coarse-grained de- 
scription plays a crucial role for the order-parameter dynamics. Coarsening is intimately 
related to the fact that domain growth is a scaling phenomenon: domain pat terns at dif- 
ferent times differ solely by a global scale factor (see the review bv lBravl2 002). A suitable 
coarse-grained description for systems where the order parameter is not conserved (e.g. for 
anti-ferromagnetic ordering) is provided by the time-dependent Ginzburg-Landau equa- 
tion. When the order parameter is conserved, as in phase sepa ration, the coarse-grained 
dynamics is ruled by the Cahn-Hilliard equations llBravll2002i) : 

^ = a 2 — (11) 

dt x 5w' [ ' 

where u>(x,i) is a suitable coarse-grained order-parameter and F a Landau free-energy 
functional: 

x \v w \ 2 + i(w) . 



F[w] = J 



dx 



_ 2 ,~, .-^j. (1.2) 

The potential I(w) typically has a double-well structure, whose minima correspond to 
two equilibrium states. A is a positive constant related to the distance between the 
equilibrium states and thus the size of the interface between them. 

In fluid mechanics, the Cahn-Hilliard equations (|1.1|) play a fundamental role in the 
stability analysis of large-scale perturbations. In a variety of situations, it turns out 
that the evolution of large-sca le perturbations is governed by equatio n (II. Ill , with a 
fourth-order potential I (w) (see iNenomnvashchvil 1 9 7fil ISivashinsk vil98,4IPedloskvll9S7t 
llvTa.nfroi fc You nd Fl 999) . The structure of the potential controls the profile and the i nter- 
actio ns of the so-called kink-antikink structures observed in snapshots of the flow l|Shel 
119871) . 

In the present paper, we focus our attention on a simple model of viscoelastic flows, 
the so-called vi scoelastic Kolmogorov flow. Its linear stability analysis has been recently 
investigated by IBoffetta et all (|2005a ), while the turbulent regime and its massive drag 



reduction effects have been studied by Bof fetta et all I 200Rbh . Here, we analyse the weakly 
nonlinear dynamics, intermediate between the linear stage of evolution and the fully 
turbulent regime. 

The starting points of our analysis are three results obtained bv lBoffetta et all (j2005a) 
for the linearized stage: i) The most unstable perturbation has a long wavelength (large- 
scale) compared to the period of the basic Kolmogorov flow; ii) Its linear evolution is 
captured by asymptotic multiscale methods, at least up to moderate elasticities of the 
flow; iii) The most unstable perturbation is transverse with respect to the basic flow. 

Multiscale asymptotic methods can be applied, as in the Newtonian case, to show 
that the evolution in the presence of polymers obeys a one-dimensional Cahn-Hilliard 
equation of the form 1)1.1(1 . The point demonstrated here is that there exists a critical 
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value of the elasticity, where the potential I(w) passes from the fourth to the sixth order 
in the field w. This corresponds to a triple critical point. Due to the one-dimensional 
character of the nonlinear dynamics, the transition strongly impinges on the dynamics 
of the large-scale perturbation. 

Above the critical elasticity, "hydrodynamic" kink-antikink structures are replaced by 
generalized kinks and anti-kinks and their annihilation processes are shown to be severely 
slowed down. Moreover, below the critical value of the elasticity, the mechanism of insta- 
bility is linear and nonlinear terms stabilize the flow. Conversely, above the critical value, 
we show that a sub-critical, nonlinear mechanism of instability takes place, provided the 
initial amplitude of the perturbation be sufficiently strong. 

The paper is organized as follows. In §0 and |2| we d escribe the viscoelasti c model 
considered in the sequel and briefly review the results by iBoffetta et al\ l|2f)05aj) needed 
here. In §0] we use multiscale methods to derive the coarse-grained equations for the 
perturbations. In §|SJ we study the system around the triple critical point and work out 
the evolution equations in its neighborhood. In §[H]and|3 we reformulate the asymptotic 
behaviour of the coarse-grained equations in terms of variational analysis and present the 
numerical results that corroborate our analytical predictions. Finally, in §[S]we address the 
problem of drag reduction and show that, even for the weakly unstable regime considered 
here, the injection of polymers induces an enhancement of the mean flow amplitude. 



2. The viscoelastic Navier Stokes equations 

Several models have been introduced (see e.g. lHincr]ll977j) to describe viscoelastic 
fluids. A powerful class describes the fluid as non-Newtonian, accounting for the reaction 
of the polymers onto the flow via an extra-term in the stress tensor. A popular and often 
employed model within this class is the Oldroyd-B llOldrovdlllflRrl . which is the one 
considered in the sequel. We briefly review it here for the sake of completeness. 

In the Oldroyd-B model, it is assumed that viscoelastic flows can be treated as a dilute 
suspension of elastic dumbbells, i.e. identical pairs of microscopic beads connected by 
Hookean springs. The flow is considered "external" to the molecule, neglecting the effects 
of the finite size of the polymers on the flow. Furthermore, the polymer concentration is 
supposed to be uniform and low enough to neglect polymer-polymer interactions. 

The reaction of the dumbbells on the fluid is treated at a mean-field level and the 
study of the dynamics is limited to scales much larger than the inter-polymer distance. 
The polymer solution is regarded as a continuous medium, whose reaction on the flow is 
described by an elastic contribution T to the total stress tensor of the fluid. Its value per 
uni t density depend s on the free energy of the molecule and the thermal noise as (see 
e.g. lBird et iHll987l) : 

T = -n p (RF) - n p k B ei, (2.1) 

where n p is the polymer density, fceO is the energy associated with thermal noise, Fi is 
the dumbbell relaxation force and Ri its elongation vector. The average is taken over the 
statistics of the thermal noise. Assuming the force between the beads to be Hookean with 
dynamical coefficient Kq, the average in the elastic stress reduces to (RF) = —Ko{RR)- 
The latter is proportional to the conformation tensor <x = (RR)/R%, where i?o denotes 
the equilibrium spring length. The inclusion of the extra elastic stress term in the Navier- 
Stokes equations leads to the following equation for the viscoelastic flow: 

d f v + (v ■ d)v = -dp + vl3d 2 v + ^—P±d ■ (a - 1) + / . (2.2) 

T 

Here, v is the total kinematic viscosity of the solution, while v/3 and /3) are the sepa- 
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rate contributions by the solvent and the polymers, respectively and we have introduced 
the dimensionless parameter (3 = rj s j '(7i p /cbOt + rj s ), t) s being the dynamic viscosity of 
the solvent, r is a parameter depending on Kq and Rq, representing the typical relax- 
ation time of the polymers. A more precise definition of r and <x comes in the following. 
Throughout the paper, it is understood that (dv) a p = d a vp and tr (dv) = d ■ v = 0. 

An equation for the dynamics of the polymer conforma tion tensor cr is needed to 
close the system of equations. Simple physical reasoning bv lBird et al\ l|l987|) gives the 
following stochastic equation for the separation R between two beads: 

R=(R-d)v-—R+\^£. (2.3) 
2t V t 

On the right-hand side, the first term is the stretching/compression term, originating 
from the spatial variation of the flow experienced by R, and the last one, is a white-in- 
time random process mimicking the effect of thermal noise on the polymers. The second is 
the relaxation due to the force between the beads, proportional to the elongation deriva- 
tive of the dumbbell free energy —dE/dRi — —d(l/2KoR 2 )/dRi. A quadratic form of 
the potential, and thus a linear Hookean force, is an approximation valid for moderate 
polymer elongations. The dynamical coefficient r is the same as the one appearing in 
l|2.2|) . Considering it constant amounts to assume that the polymers have only one relax- 
ation time. Numerical an d theoretical studie s point out that this hypothesis is reasonable 
llGeraschenko et alEoO^ . Experiments (see lLumlevlll969t IVirlJl975t iNadolink fc Haighl 
Il995h show that polymers have a spectrum of typical relaxation times, but they also 
show that interactions with the fluid mostly depend on the largest one, that is the one 
we are retaining. 

Multiplying lf!H|) by R and averaging over the statistics of the thermal noise the 
following evolut ion equation for the conformation tensor a = (RR) / R 2 is obtained (see 
iBird et »/.lfl987li: 

d t (T + (v ■ 8)a = (dv) T (t + ct- (dv) - -(tr - 1) . (2.4) 

T 

Summarizing, the set of equations (|2.2(l and (|2.4() constitutes the Oldroyd-B model 
that we shall be considering in the sequel. 

Our first step in the investigation of the effect of polymers onto the stability of the 
flow will be to find out the basic equilibrium state. The state will then be perturbed and 
the resulting equations analysed using multiscale methods. 

2.1. A basic equilibrium state 

Finding analytically the basic equilibrium state for a generic forcing / is a hopeless 
task already for the Navier-Stokes equations without polymers. The problem is further 
complicated here by the additional term in (|2.2|l and the coupling with (|2.4H . 

The problem simplifies for / = (f(z),0, 0), inducing a parallel flow U = (U(z),0, 0), 
which trivially annihilates the advective nonlinear term in (|2.2|) . A further substantial 
simplification comes from the viscoelastic version of Squire's theorem (see Appendix |A"|) . 
stating that, for parallel flows, the most unstable perturbations are two-dimensional. We 
shall the refore restrict to a two -dimensional flow (u x , u z ), without any lack of generality 
(see also iBoffetta et qi.ll2005a). We further assume f(z) = F cos(z/L), producing the 



well-known Kolmogorov flow l|Arnold & Meshalkm1ll96rf) U(z) = (V cos(z/L), 0), where 



V = Fq L 2 1 v. The corresponding conformation tensor at equilibrium is: 

(l + 2r 2 (d z Uf rd z U \_ ( 1 + 2r 2 £ sin 2 (f ) -r£sin(f)\ , 
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This choice also allows to precisely define the Reynolds number of the flow as Re = VL/v. 
In this model the elasticity of the polymers is taken into account by the relaxation 
time r only. We can thus introduce an adimensional parameter, the Deborah number 
De = tV/L, to characterize the elastic properties of the flow. 



3. Some previous results on the linear stability analysis 

It has long been known t hat the Newtonian Kolm ogorov flow becomes unstable for 
Reynolds numbers Re > \J~2 ijMeshalkin fc Sinallll96lj) : the evolution of large-scale per- 
turbations is formally described by an effective diffusive dynamics and instabilities are 
associated to the loss of positive-definiteness of the eddy- viscosi ty tensor. 

In the presence of polymers, performing a multiscale analysis ((Bensoussan et q,Llll97St 
IBavlv et aZ.lll98^) on the linearized Oldroyd-B model, one obtain s an explicit expressio n 
for the eddy- viscosity tensor, valid for sufficiently low elasticity (Bo ffetta et aLll2005a[) . 
The resulting stability curve in terms of the Reynolds and the Deborah number is re- 
ported in figure n 

3.1. Multiscale analysis 
Substituting v = u + w into H2.2I2.4J1 . the equations for the perturbation fields read: 

dw = 0, (3.1) 

dtw + d ■ (uw + wu + ww) = —dq + vf3d 2 w + v (1 — j3) r" 1 d ■ C ■ (3-2) 

d t C + d- {uC + wcr + wC) = {duf ■ C + {dwf ■ cr + (dw) T ■ C + 

+C ■ {du) + cr ■ {dw) + C ■ (dw) - t^C , (3.3) 

where q and £ are the perturbations associated to the pressure term p and the basic stress 
tensor cr. In the linear stability analysis, the nonlinear t erms containing the p roduct of 
two perturbation fields are supposed to be negligible (see iBoffetta et aLll2005a|) . 

As in the Newtonian case, it is assumed that the first unstable perturbations have 
periodicity much larger than th at of the basic flow. T he validity of this assumption 
has already been investigated bv lBoffetta et all l)2005al) and is satisfied in the range of 
parameters considered here. 

In addition to the usual "fast" space/time variables x,t, describing the basic flow, 
multiscale techniques introduce "slow" variables x = ex, t = e 2 t, to describe the large- 
scale flow, and prescribe to treat the two sets as independent. This leads to the expansion 
of the differential operators: 

d t - d t + ed h d t ^d t + e 2 d t , (3.4) 

and of the fields: 

w = (z, t, x, z, t) + en/ 1 ) (z, t, x, z, t) + e 2 «/ 2 ) (z, t, x , z , t) + . . . , 

q = q(.o) ( Zj £ ; x, z, t) + eg' 1 ' (z, t, x, z, t) + e 2 q^ (z, t,x,z,t) + . . . , (3.5) 

C = C (0) 0, t, x, z, t) + eC (1) (z, t, x, z, t) + e 2 C (2) (z, t, x, z, t) + . . . . 

All of the functions have the periodicity of the basic flow and are independent of x. 

Inserting l|3.5|l into H3.1JI - H3.3J) an d collecting terms of the same order in e, the coarse- 
grained equation for a large-scale perturbation is obtained as the solvability condition at 
the order e 2 . In terms of the stream function the perturbation evolves according to 
the non-isotropic diffusion equation: 

d t A* = v a pd 2 J}y, (3.6) 
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Figure 1. The linear stability diagram for j3 = 0.769. Stable and unstable regions are denoted by 
S and U, respectively. The bullets represent direct numerical simulations (DNS) of the complete 
Oldroyd-B system of equations, confirming theoretical predictions for this window of parameters. 

where the eddy-viscosity tensor v a p is not positive-definite for Re > \2 (in the absence 
of polymers). In general, there exist critical values of Re and De where perturbations 
start growing in time. 

The phase-space {Re, De) is thus divided in regions where the eddy-viscosity tensor is 
positive-definite (the system is linearly stable with respect to any small perturbation) 
and where there exists at least one unstable mode, as shown in figure for low Deborah 
numbers. The diagram reveals two kinds of instabilities. When the Deborah number is 
sufficiently low, the flow experiences hydrodynamic-like large-scale transverse instabili- 
ties, captured by multiscale analysis. In this region, the critical Reynolds number where 
the flow beco mes unstable, g r ows wi th De: polymers stabilize the flow. This has been 
interpreted by iBoffetta et al\ l|2005a|) , and will be shown in §[8] to be a prelude to the 
drag reduction effect observed in the turbulent regime. 

For high values of the Deborah number (not shown in figure ^) the multiscale anal- 
ysis predicts the flow to be unstable, even for very low Reynolds numbers. However, 
numerical simulations show that the assumption of scale separation does not hold and 
multiscale techniques are not applicable. This region, possibly characterized by purely 
elastic instabilities, will not be the concern of the present investigation which focuses on 
< De < 2. 

If the amplitude of the large-scale perturbation is strong enough and/or the eddy- 
viscosity is negative, nonlinear effects are important and should be taken into account. 
These two situations correspond to different scalings of the fields and will be treated 
separately in the next sections. 

4. Nonlinear dynamics: the standard Cahn Hilliard equation 

Linear stability analyses, by their very definition, are not able to capture the full-time 
dynamics of unstable perturbations: as perturbations grow in time, their magnitude 
becomes large and nonlinearities ought to be taken into account. In the Newtonian case 
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l|Sivashinskvlll98.4 IShelll987h . the lowest-order nonlinearity (third-order) is sufficient to 
stabilize the linear exponential growth and lead to a steady state. Here, we show how 
this occurs and generalize it first to the viscoelastic case. The next section will then be 
devoted to the case where the lowest-order nonlinearity does not stabilize the flow and 
higher-order nonlinearities become relevant. 

A careful analysis of the linear eddy-viscosity tensor v a p derived in (|3.6[) ensures that 
for low enough Deborah numbers the first modes to become linearly unstable are the 
large-scale transverse ones. For barely unstable flows we may expect that the perturba- 
tions involved in the nonlinear dynamics will be confined to these modes. This suggests 
that the result will be a one-dimensional diffusion equation for the averaged transverse 
modes, linearly stable for small-scale modes, involving at least one nonlinear term. 

Assume now that the initial amplitude of the large-scale perturbation is sufficiently 
small and the system we consider is in the surroundings of a point of the critical curve. 
According to l|3.6|) . the average transverse velocity perturbation (w z ) linearly evolves 
according to a diffusion equation: 



dt 



-Ad 2 



(4.1) 



where A is a positive coefficient representing the linear eddy- viscosity tensor of (|3.{jl) . 
restricted to low Deborah numbers (as explained in the end of §[3J). ft vanishes on the 
stability curve, being positive above it and negative below it. 

In this equation, all modes are linearly unstable. It needs to be modified to keep track 
of the multiscale hypothesis, which requires small-scale modes to be stable. This is done 
introducing a bi-Laplacian term into l|4.2(l to stabilize the small scales (the fourth-order 
derivative ensures that this term will be dominant on the small-scale perturbations only) : 



-Ad 2 



ca 4 « 



(4.2) 



In general, close to the linear instability threshold, where the coefficient A changes sign, 
we do not expect C to vanish. To comply with the stability requirements, we will request 
it to have a finite, positive value in the region of interest. 

We expect to find the presence of a nonlinear term, eventually stabilizing this growth. 
This part cannot be played by the advective nonlinearity because of the one-dimensional 
character of the equation. The next-order nonlinearity is cubic and must contain at least 
two space derivatives: one before the whole term, to ensure momentum conservation, 
and an additional one to respect space-inversion symmetries. This yields a nonlinear 

term: Bd x (^{w z 1 ^) 2 d x ('w z 1 ^)^j , where B is some constant related to the (nonlinear) eddy- 
viscosity. 

We can now introduce, as in §03 the "slow" variables x = ex and t (notice that we still 
do not know the scaling between t and t). The space derivatives must again be expanded 
as di — + di + ed\ . We now have to look for a prescription on how to expand the different 
fields in terms of e. 

In the vicinity of the marginal eddy-viscosity curve A ss 0, and a Taylor expansion 
gives v — v c ), where v c indicates the critical viscosity. Balances between the 

term Ad(w z ) and both the cubic nonlinearity and Cd 4 (w z ) yield: 



Be 2 el 



dA 



)e 2 e„ 



dA 

dv 



{v - v c )e 2 e v 



Ce 4 e z 



(4.3) 



Here, e w is the scaling of the amplitude of the velocity perturbation w. 

Equations Ij4.3|l completely define the scaling for the velocity perturbation and the 
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distance from the critical viscosity: 

e = e w , »—^~e^ V = Vc {l-e 2 ). (4.4) 

It follows from 14.41) that the Reynolds number Re — Re c (l + e 2 ). The comparison of 
any of the previous terms with the (slow) time derivative of w gives the scaling t = e 4 i. 
As for the scaling of the polymer conformation tensor, balancing £/r and (dw z ) ■ a, we 
obtain that the scaling of C coincides with e w . The same equality holds for the pressure 
field. 

Summarizing, the fields are expanded as: 

w = ew^ (z, x, t) + e 2 u/ 2 ) (z,x,t) + . . . , 

q = eq<H (z, x, i) + eV 2) (z, x, t) + . . . , (4.5) 
C = eC (1) (z, x, t) + e 2 C (2) (z, x, t) + . . . . 

The next step to obtain a coarse-grained equation for the large-scale dynamics is 
to plug (|4.5() into H3.lfH3.3fl . Exploiting the chain rule, the definitions of x and t and 
averaging along z, we end up with a set of equations involving solely the large-scale fields. 
The equation for the large-scale transverse perturbation (w^)(x,t) is obtained from 
the solvability condition at order e 5 . For details on the Newtonian case and solvability 
conditions, see iGama et al\ l|l994[) - 

We can summarize the whole procedure in the following schematic way: 

(a) Solve the continuity equation. The explicit expression of u4™ is thus obtained in 
terms of known functions of z. 

(b) Solve the equation for Q™ this can always be done algebraically as is slaved 
to the field. 

(c) Solve the evolution equation for w z . This field is obtained from (a); we are then 
able to obtain the expression for the pressure field perturbation q. 

(d) Solve the system for C^z and u4™^ by direct integration. 

(e) Algebraically obtain the explicit expression for . 

(f) Impose solvability condition at order n + 1 on the continuity and the velocity field 
equations. Such condition is automatically fulfilled by the polymer conformation tensor, 
as it is slaved to the velocity field at the previous order. 

The final equation has the form of a "standard" Cahn-Hilliard equation: 

dM 1) )=d x [(-A + B( w ^f)dM 1) )}~Cd i x (w^) . (4.6) 

"Standard" is meant to stress that the structure of l|4.6fl (i ncluding the cubic non- 
linearity) emerges in a variety of hydrodynamic situation s (see iNepomnvashchvil 1 1 9 76t 
ISivashinskvlll985l |Pedloskvlll987t iManfroi & Youndll999|) . The parameters A,B,C are 
known functions of the parameters De and /3, as shown in figure |3 It is worth noting 
that A is non-negative as the system is supposed to be slightly above the threshold of 
instability and we have explicitly incorporated a negative sign in (|4.6fi . 

The saturation of the instability requires two conditions. First, C must be positive 
to ensure that the instability be saturated at sufficiently high wave-numbers (still much 
smaller than those of the basic flow, of order unity). Second, B ought to be positive to 
ensure that, as {wi 1 ^) becomes 0(y/ A/B), the nonlinear eddy-viscosity —A + B{w z 1 " 1 )' 2 
change sign and the growth be again saturated. Both these conditions are satisfied up to 
a critical value of the Deborah number, De* (see figure EJ). 
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Figure 2. The parameters A,B and C appearing in the coarse-grained Cahn-Hilliard equation 
114.611 as a function of the Deborah number (for /3 — 0.769). 

Up to the critical Deborah number De* , the equation (|4.6|l for the large-scale pertur- 
bations has the same structure as in the Newtonian case and the only difference is in 
the numerical value of the parameters A, B, C . As we shall see in the next sections, this 
property ceases to be true above De* . 

To conclude, we stress the fact that all the fields up to order four are expressed in 
terms of explicit functions of the fast variables and of the large-scale field {w^P), obeying 
the Cahn-Hilliard equation 1)4. 



5. Generalized Cahn Hilliard dynamics 

We have observed in the previous section that, along the marginal linear stability curve, 
there exists a critical value of the Deborah number, De* , where the cubic nonlinear term 
becomes negative. Furthermore, the change of sign is taking place in the region where 
the small-scale operator is stable. The problem is thus well-posed and lends to multiscale 
methods. 

For De > De*, the field keeps growing at sufficiently large scales, until it reaches ampli- 
tudes where the next-order nonlinearity becomes important. Its structure is dictated by 
the conservation of momentum and the symmetries of the basic flow: d x (^(w z ) 4 B x (w z )^j , 
with a regular coefficient D in the neighborhood of the critical point P* , where both the 
eddy-viscosity and the coefficient of the third-order nonlinearity change sign. 

Four regions can be identified around P* (see figure[3J|. The eddy- viscosity A — curve 
has been obtained by means of the linear stability analysis (§[3J) . The linear approximation 
of the curve B = in the vicinity of P* is obtained from the analytic expression of B on 
the marginal curve and the marginal curve itself. 

Zone I is linearly unstable (A > 0), has a third-order destabilizing term (B < 0) and 
we can guess that a fifth-order term will enter into play to stabilize the growth. Zone II is 
particularly interesting as it is linearly stable (A < 0), but has a third-order destabilizing 
contribution (B < Q). Perturbing with a field of sufficiently strong amplitude, the system 
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Figure 3. The phase-space around the critical point P* where both the eddy- viscosity and 
the coefficient of the third-order nonlinearity change sign. The region is divided in four regions 
schematically sketched here by the two critical curves A = and B — 0. The former is found 
from the linear stability analysis in section 3. The latter is found locally, around the ^4 = 
curve, by solving l|5.5^ . and is linearly extrapolated for graphical purposes as a dashed line. For 
P = 0.769, the curve B — is inclined at approximately 60° with respect to the De axis. 

jumps to the asymptotic steady state where the two nonlinear terms (third and fifth- 
order) balance each other. Zone III is completely stable (A < 0,B > 0). In the last 
region, IV, as De approaches the critical value, the coefficient B goes to zero and cannot 
saturate the exponential growth from the linear instability. The fifth-order nonlinearity, 
which is negligible far from the critical point, must enter again into play. 

5.1. Zone I 

When both the Reynolds and the Deborah numbers exceed their critical values, previous 
considerations suggest the following structure for the coarse-grained equation: 

d t w = -Adlw ~ \B\d x {w 2 d x w) - Cd^w + Dd x (w 4 d x w) . (5.1) 

Confining the analysis to the surroundings of the critical point P* , we may represent 
the position in phase space as: 

v=v*(l-K 1 e l/ -K 2 el), (5.2) 
De= De*(l + e Dc ). (5.3) 

Adequately choosing the K\ and parameters, any point around P* can be reached 
as e varies. The reason why we need to incorporate in (|5.2J) the additional contribution 
of order e 2 will be clear shortly. 

In the neighborhood of P* , the coefficients A and B are expanded as: 



OA 



dA 



B = 



dDe 
dB 

dDe 



{De-De* 



dv 
dB 

dv 



(5.4) 
(5.5) 



where all derivatives are computed at P*. 
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The scaling in e of the velocity field amplitude, e w , and the parameters e„, e De is found 
by requiring that all terms in l|5.1(l be of the same order in the scale-separation small 
parameter e. 

The comparison between the last two terms in 1)5. 1[) fixes the relation between e and 



D e 2 ei ~ C e 4 e w e v 



,1/2 



(5.6) 



The parameters e„ and e De are found by comparing the terms associated to A, B and 
D in |JOJ. Using we obtain: 



oL>e cw 
O £ 2 e 5 / 2 ~ [^-(e^e*) - e 2 e 3 / 2 . (5.8) 

Choosing e„ = e De = e and setting ifi to ensure ~§]j{>De* — ^K\V* = 0, both 

equations 1)5. 7|) and (|5.8|) are satisfied. Equation (|5.7II is balanced by the second-order 
term of the v expansion l)5.2p . dependent on if 2. The scalings of time, pressure and 
polymer conformation tensor perturbation, e 4 , e 1 / 2 and e 1 / 2 , respectively, are derived as 
in §H 

Once the scalings have been determined we can proceed as in §01 to obtain the large- 
scale equation for (wi 1 ^ 2 ^)^, x). The evolution equation emerges now from the solvability 
condition at the order e 9 / 2 : 

d t (w{ 1/2) ) = d x [{-A + B{w^f+D{w^t)d x {w^^ (5.9) 
where the coefficients are explicit functions of (3. For f3 = 0.769, they read: 



A = 0.5106 + 1.965^2 , B = -8.979, 
C = 0.9439, D = 23.11, K ± = 0.594. 



(5.10) 



Although (|5.9|) belongs to the class of the Cahn-Hilliard equations l|l.l|l , the emergence 
of the new, sixth-order nonlinearity will be responsible for new dynamical aspects, not 
present for De < De* , which will be discussed in detail in §|SJ 

5.2. Zone II 

For Deborah numbers above the critical value, perturbations are nonlinearly unstable: 
B < 0. This is true regardless of the sign of the linear term and strong enough pertur- 
bations may then grow even if the system is linearly stable. 

Let us then consider systems with v > v* and De > De* . No major difference with 
respect to case I is expected. At zero-th order, the coefficients A and B vanish and equa- 
tions l|5.4|) - (|5.5|) hold. Again, we define the position in phase-space via the two parameters 
e„ and ejj e - As the viscosity is now larger than the critical value, a positive sign appears 
in the expansion of the viscosity: 



(5.11) 



while l|5.3|l holds. The parameter K2, as we shall point out later, can take any value 
compatible with the condition A > 0. 

The same calculations discussed in the previous subsection can be carried out to derive 
the coarse-grained equation for the transverse velocity. As one might expect, its form is 
exactly the same as (|5.9|l . a generalized Cahn-Hilliard equation. The only difference is 
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in the value of the parameters. For j3 = 0.769, they read: 



A = -0.2202 
C = 0.9439, 



- 1.965^2 , 
D = 23.11 



B = -35.62, 
K x = 0.5974 . 



(5.12) 



Only A and B have changed with respect to l|5.10|l . as expected since they are the only 
parameters which depend on e (and thus on Re and De) in physical coordinates. Notice 
that there is an upper bound on the values we can choose for reflecting the linear 
stability requirement. 



5.3. Zone IV 

What happens when the Deborah number is barely smaller than the critical value De*l 
Sufficiently close to it, the third-order instability can be made subdominant with respect 
to the fifth-order and our aim here is to work out the scaling coefficients corresponding 
to such situation. 

For this purpose, let us assume that the cubic nonlinearity is negligible. At leading 
order, the terms associated to A, C and D must be of the same order. This means: 



4 

e e„ 



e 2 e 5 



e 2 e„ 



e 2 e,, 



OA 
We 
' dA 



dDe 



(De-Be*) 
(De-De*) 



dA, ' 
av 



(5.13) 
(5.14) 



and implies: 



v = v* (1 - K 2 e 2 ) , De= De*(l-e 2 ). (5.15) 

Additionally, the velocity field scales as e 1//2 , as the pressure and polymer fields do. The 
time derivative scales as e 4 . 

To be consistent, we are left to check that the third-order nonlinearity is negligible. 
Using the previous scalings and the ensuing fact that B ~ 0(e 2 ), we have to verify that: 

0(Bd 2 w 3 ) < 0(Dd 2 w 5 ) 0{e 11/2 ) < 0(e 9 / 2 ) , (5.16) 

which holds true. It is now possible to apply the strategy discussed in §0]to derive the 
large-scale equation and obtain (at order e 5 ): 



~d t {w£ /2) ) = 9 X [(-A + Diwi 1 ^) d x (w£ /2) ) 



Cdiiw 



(l/2)\ 



(5.17) 



where C and D have the same value as before and A = 1.1740 + 1.965^2. 



6. Variational formulation 

It is well-known that the "standard" Cahn-Hilliard equation admits a var iational 
formulation in terms of a Ginzburg-Landau potential iCahn fc rlilliardlll95£ft . Equa- 
tion (|4.6() . after appropriate rescalings, w — > (A/BY^w, t — > A~ 1 t,\ = C/A, is recast 
in the form l|l.l|) with the Lyapunov functional: 



F\w] = 



±(d x w) 2 +I(w) 



2 4 

1 Tt \ W W 

dx, I(w) = -- + -. 



(6.1) 



Note that mean fields only are considered, that is w must be read as the rescaled leading 
contribution (wi 1 " > )(x, t). 

The existence of a Lyapunov functional implies the existence of an asymptotic state 
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for w, if boundary conditions are periodic and stationary. Such state corresponds to a 
minimum of the Lyapunov functional and it is calculated by the following equations: 

I'(w) = Xd x w <-> d x I = ^d x (d x w) 2 . (6.2) 

When w is a maximum (or a minumum), I(w max ) is constant and w max is obtained 
solving I' = as w max = ±\/3- Considering this boundary condition, equation l|6.2|) can 
be easily solved. Its solutions are the well-known kink and anti-kink structures, namely: 



10 



±V3 tanh 



2X X 



(6.3) 



The issue now is whether or not a Lyapunov extremal formulation exists in the gener- 
alized Cahn-Hilliard case 15.111 a s well, and how i t relates to the standard one. In partic- 
ular, a Painleve test llAblowitz fc Clarksonlll99ll) can be performed on the equation to 
check its integrability. The calculation consists in checking that all movable singularities 
(whose location depends on initial and/or boundary conditions) are poles (see for details 
lAblowitz fc Clarksonlll99l . The test is based on a well-known connection between the 
integrability property of a nonlinear differe ntial equation and its analytic structu re for 
complex values of the independent variable l|Kowalesvkilll889l Il890t rPainlevelll897|) . The 
explicit calculation is performed in Appendix^] The generalized Cahn-Hilliard equation 
indeed enjoys the Painleve property and is thus integrable. 

Let us then write the equation (|5.1I) after the rescalings w — > (A/ BY^w, t — ► 
A~H,X = C/A,<y = AD/B 2 : 



d t w = -flgui - \d 2 w 3 - Xdjw + ^dlw 5 . (6.4) 

Integrability of this equation is related to the existence of the following Lyapunov func- 
tional, similar to that of the standard case, yet with a sixth-order potential: 



F\w] = 



~(d x w) 2 + I(w) 



dx , ICw) = — + -^w 6 . (6.5) 

v J 2 12 30 v ; 



The calculation of the function corresponding to its minimum is performed using again 

era . 

All solutions tend to final steady states which minimize F. The approach to the so- 
lution is however nontrivial and the structure is made of plateaux having velocity ±Wq 
(/'(Wo) = 0), separated by positive and negative kinks (see figure^}. The amplitude of 
the velocity w in the plateaux is: 



5 + 725 + 180^ 

67 



Note that, at small 7's, the asymptotic velocity W diverges as 1/^/7. This is quite 
intuitive: the field amplitude equilibrating the third and the fifth-order nonlinearities 
increases as the coefficient of the fifth-order nonlinearity reduces. 

The explicit expression of the profiles for kinks and anti-kinks is obtained from the 
integration of equations l|6.2|l . (|6.5|l . For the sake of example, when A = 1/2 and 7 = 10/9 
the profiles read: 

2</3x _ ^ 

w = ±V15 - . (6.7) 

V5e 4 ^ + 26e 2 ^ + 5 
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Figure 4. The "generalized" (solid) and "standard" (dotted) kinks for A = 1/2 and 7 = 10/9. 
The former has a manifestly shorter range. It is shown in the body of the text that this entails 
longer time-scales for their annihilation with the corresponding anti-kinks. 




Figure 5. The potentials associated to the different evolution equation. Curve a is related to 
the standard Cahn-Hilliard equation (fourth-order potential) ; curve b represents the generalized 
Cahn-Hilliard equation (sixth-order potential) . Curve c is the characteristic triple- well potential 
of the purely nonlinearly unstable case. The plots are in arbitrary units, to ease the comparison 
between the curves. 



The generalized kink-antikink structures, e.g. those given by equation l|tj.7|) . will be 
dubbed "generalized" kinks and anti-kinks. 
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6.1. Dynamics of generalized kink/antikink annihilation and approach to equilibrium 

Detailed calculations are performe d following lLegras fc Villond l)2003|) . who in turn based 
theirs on iKawasaki fc Ohtal l|l982|) . They are lengthy, yet quite simple in their basic idea, 
sketched hereafter. 

During metastable transitions, the kinks do not satisfy (|6.2[) exactly, due to the pres- 
ence of other kinks and/or anti-kinks. The deviation of the amplitude in the plateau is 
proportional to e~ sA , where A = 4\x\ and x denotes the distance to the point where 
w = 0. Here, s is the inverse of the typ ical length scale of this deviation (for details, see 
Appendix A of lLeeras fc Villonejl2003|) . The quantity s turns out to be crucial as neigh- 
boring kinks and anti-kinks attract proportionally to e~ sAx , where Ax is the distance be- 
tween neighbouring kinks and anti-kinks (for details, see Appendix B of lLeeras fc Villonel 
120031) . 

The behavior of the kink size s is grasped as follows. Consider a metastable state of the 
Cahn-Hilliard equation. The potential felt by a kink w(x) close to the plateau w = Wo 
is estimated by the Taylor expansion: 

I(w - W ) ~ I(Wo) + I'(W )(w - Wo) + I"(W ) {W ~^° )2 , (6.8) 

where we know that I' (Wo) = 0. Note also that the dynamics of w does not change if we 
add an arbitrary constant to the potential /, so that we can set /(Wo) = 0. 

Let us now calculate the shape of the profile between w and Wo- For a metastable 
state, d t (w — Wo) = 0, that implies: 

^(d x (w - Wo)) 2 + [I"{W ) {W ~Y° )2 ] = • ( 6 - 9 ) 

Interpreting d x as the inverse of the typical length scale s for (w — Wo), one easily obtains 
s = \/ 1" (Wo). The second-order derivative can be explicitly calculated using Ij6.6|) : 

I" (Wo) = 4+^W 2 . (6.10) 

Qualitative properties of s are easy to grasp. At large 7's, the size of the kinks tends 
to a constant, independent of 7. At small 7's, the kinks get steeper and steeper, their 
size scaling as 7 1 / 2 . This implies that the convergence to equilibrium will be slower and 
slower as 7 is reduced (recall that the kinks attract proportionally to e~ sAx ). 

For the same band of unstable modes, i.e. keeping A fixed, it holds that the convergence 
to equilibrium is slower for the generalized than for the standard Cahn-Hilliard equa- 
tion. Indeed, for the Cahn Hilliard potential Ich = —w 2 /2 + w 4 /12, the second-order 
derivative I'c H (W ) = 2. As for (|jHU|) . we can use the identity 1 + W 2 /3 = 7W 4 /5, 
following from the very definition I'(Wo) = 0, to obtain /"(Wo) > 2. This implies that 
the interactions for the generalized kink-antikink structures have a shorter range and 
their dynamics of annihilation is thus slower. 

A special remark applies to the linearly stable case (zone II). In this case, the equation 
is associated to an uncommon triple-well potential. The typical nonlinear kink-antikink 
dynamics appears only if the initial perturbation will be energetic enough to let the 
system "jump" out of the central well and fall into one of the side wells. 



7. Numerical results 

The analytical results presented in this work have been obtained by multiscale tech- 
niques. Their basic assumption is the strong scale separation between the basic flow and 
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Figure 6. Growth rates g of the transverse Fourier modes k for a simulation with De = 1.4 and 
(3 — 0.769. The simulations are performed in a rectangular domain with aspect ratio 1/64. The 
distance to the critical point Re/ Re c — 1 is 0.28. The solid line represents the linear prediction 
17. II . The circles representing the numerically computed growth rates have been obtained with 
a DNS simulation by a linear fit of the logarithm of the energy for each mode versus time, in 
the early stages of their exponential growth. 

the most unstable perturbations. In this section, we shall present the numerical simula- 
tions performed to check the validity of this assumption. 

The linear analysis results have already been checked in lBoffetta et all ()2005a|) by re- 
ducing the original linear partial differential equation to a generalized eigenvalue problem 
and computing its eigenvalues/eigenvectors. The scale separation hypothesis is found to 
be well verified up to Deborah numbers of order unity (De « 2.3 for the value (3 = 0.769 
used in this study). For larger De, scale separation does not hold and multiscale methods 
are not applicable. 

To check the nonlinear results derived here, we have numerically simulate d the complete 
Oldroyd-B model equations H2.2J) - H2.4[) via a pseudospectral algorithm fsee lCanuto et all 
1988, for details on the numerical method). In the following, we will refer to these as to 
Direct Numerical Simulations (DNS), while numerical integration of the one-dimensional 
Cahn-Hilliard equation (|4.6I) will be referred to as CH simulations. 

To enforce a transverse perturbation, we integrated the equations on a rectangular 
slab with L = L x = 2ir and L z = 64L X . The aspect ratio r is then fixed at 1/64. The 
simplest check of our results concerns the growth rates of the instability which, in the 
linear regime, can be obtained by the Cahn-Hilliard equation. Neglecting the nonlinear 
term, the dispersion relation for the transverse Fourier modes k reads: 

Re - * -.2 r.L4 



In figure we report the growth rates of the first modes for a (white-noise in space) 
small initial perturbation. We are then able to observe also negative g (stable modes). The 
comparison with the linear prediction is excellent, even for modes whose scale separation 
is not very small. 
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Figure 7. The energy associated to the lowest wavenumber modes resulting from a CH sim- 
ulation. The values of the coefficients have been arbitrarily chosen for conveniency of display. 
The quasi-stationary states can be clearly seen up to the asymptotic one corresponding to the 
largest periodicity. In this simulation, De = 1.4 and Re/Re c — 1 = 0.28. 

Let us now consider the nonlinear stage of the perturbation growth. It is well known 
that the time evolution of the Cahn-Hilliard equation shows a succession of long-lasting 
metastable states characterized by a well defined periodicity. For sufficiently small initial 
perturbations, the wave- number k associated to the maximum growth-rate g will be the 
first to reach the balance between the destabilizing linear term Ad 2 (w z ) and the stabiliz- 
ing non- linear one Bd 2 (w z } 3 . When such equilibrium is reached, the energy associated to 
that mode is constant and the system is quasi-stable. In the meanwhile the other modes 
k m ax — 1, k max — 2, . . . keep growing. When the mode k max — 1 balances the two terms, 
the energy associated to the mode k max drops. This new state is again quasi-stationary 
and has a well-defined periodicity k max — 1. 

The process continues until a state with the box periodicity is reached (see figure [£} ; 
such a state is stationary and corresponds to the asymptotic behaviour in §0 The kink 
structures described there, are characteristic of all of these stages. Indeed, any transition 
between two quasi-stationary states can be seen a s a kink-antikink annihilation, yielding 
a decrease in periodicity, as in figure 181 (S hdll987|) . 

To check the results obtained in §0] we have performed a DNS simulation for a partic- 
ularly long lapse of time. The excellent agreement between the DNS and the prediction 
of the Cahn-Hilliard equation is shown in figure [UJ 

The same comparison can be realized in the neighborhood of the critical point P* . 
This kind of simulation is much harder than for the standard Cahn-Hilliard, because it 
involves a very precise knowledge of the position of the critical point, and there is no easy 
way to obtain this from the simulations. Moreover, any system we simulate will be at a 
finite distance from the critical point. The parameter that will mostly feel this difference 
will be D, as we have chosen it to be approximately constant around P*. We have been 
able to overcome this weakness via a limited tweaking of the D parameter in the CH 
simulation. As shown in figure 1101 an excellent agreement between the curves is again 
achieved. 
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Figure 8. Instantaneous transverse velocity field at different times. The simulation is the same 
as in figure The transition between two metastable states can be regarded as a kink-antikink 
annihilation. In this figure a transition from afc = 2toafc = l state is represented, the a;-axis 
being the physical x direction of the integration box and the y-axis being the amplitude of the 
w perturbation. The time figure set over the graphs refers to the evolution shown in figure |7| 
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Figure 9. The comparison between DNS simulations of the Oldroyd-B model and the 
coarse-grained Cahn-Hilliard equation derived in the body of the text. The thicker lines rep- 
resent the evolution of the lowest-energy modes in a DNS simulation, while the thinner lines 
are the result of a CH simulation. Its dynamical parameters have been set with the results ob- 
tained in section [I] This particular figure refers to a simulation with De = 1.4, /3 = 0.769 and 
Re/Re c - 1 = 0.14. 



8. Clues on drag reduction 

One of the most striking properties of viscoleastic fluids is the drag reduction effect. In 
1949, Toms found that the injection of minute amounts of polymers in turbulent fluids 
flowing in a channel was able to increase the mean flow by an amount soar i ng up to 80% . 
Even if this phenomenon has been known for over fifty years llTomslll949t iLumlevlIT^ 
a satisfactory understanding of its fundamental mechanisms is still lac king. 

A l arge number of experiments has been performed to s tudy this effect (see, e.g. IVirkl 
ImaiNadolink fc Haighlll99.'Tt ls reenivasan fc Whitel fcoOCn. but a burst in its theoretical 
analysi s occured after drag redu ction was found in numerical simulations of viscoelastic 
fluids l|Sureshkumar et alJ ll997j) . T he activity i s being spurred both by fundamental 
interest and industrial applications l)Larsonlll992(K 

Drag reduction is commonly associated to channel flows and boundary effects. Still, it 




Figure 10. The generalized Cahn-Hilliard equation reproduces the dynamics of the Oldroyd-B 
model around the critical point P* . As in figure|5| the thicker lines are DNS simulations while the 
thinner ones are CH simulations. This comparison was realized for De = 1.62 and Re — 2.5159. 

is now clear that the phenomenon is present even for free flows l|Boffetta et o,/j l2005b'l. 
What we show here is that, even at relatively small Reynolds numbers, an increase in 
the Deborah number produces an enhancement in the mean flow amplitude. Simply 
looking at the linear stability diagram JQ) we may already conclude that, as the polymer 
elasticity grows, so does the critical Reynolds number and the flow is stabilized. Let us 
further investigate this effect analytically using the results of §0] 

A parameter that c an be employed to stu dy the mean flow properties in free flows is 
the drag coefficient / jBoffetta et gi.ll2005bh : 

The drag coefficient can be seen as the ratio between the energy input (through the 
forcing Fq) and the mean energy of the flow. As we are interested in mean effects only, 
we will average U 2 over the basic flow periodicity. This will ensure that only mean effects 
will be taken into account. 

When the state is linearly stable (low Reynolds numbers) we know that no perturbation 
can alter the basic flow, U = V = FqL 2 /v and thus / = Re~ x . 

In §0J we have solved all the equations of motion up to the fourth order. They give 
the following form of the flow (up to the second order): 



V(L 2 + {fi - 1)ut) 



z . 



U x {z) = Vcos(-) + ^ ' ^ ~'"' (w{V) sin (~) + (8.2) 

_^ + ^-l )( 2^ + ^)] , 
v l rL L 
DeL{B-l),~ . m ,, . ,2z N 

The first term is the basic, stationary Kolmogorov flow. Averaging over all possible initial 
conditions, component proportional to sin(z/L) and sin(2z/L) disappear. The resulting 
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Figure 11. The drag coefficient / plotted versus the Deborah number De at constant Re. As 
the polymer elasticity grows, the drag coefficient diminishes. This implies that the mean flow 
grows with De. 



expression for the mean flow in the x direction reads: 



U x (z) = {V + h{De,/3)^ F ^)cos(-) = V ef fCos{-), (8.3) 



where the quantity (wi^) 2 follows from the Cahn-Hilliard equation in the stationary 
state: 



= -8*{wP)Ae 2 + ^dHw^f - di(w^)C. (8.4) 



As (wi 1 ^) is periodic, we can integrate twice over the domain and notice that, on the 
plateau, the last term is zero. The field amplitude must then satisfy: 



o B (i) (x) 3e 2 A 



= -Ae 2 + -(w^) 2 (w^) = y ^- (8.5) 

Since the analytical expression of A and B is known, as well as how e changes with De 
for a fixed Reynolds number, the analytical expression for / is obtained: 

f= vFL = V 2 = 1 , s 

V 2 ff ReV 2 ff Re{l + h ^Re^ ? > 

where h, A, B and Re c are explicit functions of the Deborah number and (3. 

As we want to investigate how the polymer elasticity affects the flow, a meaningful 
approach is to keep the Reynolds number fixed, while varying the Deborah number. This 
allows studying how the same flow reacts when different kinds of polymers are injected. 
Once p and Re are chosen, it is possible to plot / versus De on the basis of analytical 
results, as in figure ITT1 The drag coefficient is clearly decreasing with the Deborah number 
even though the flow is barely in its nonlinear regime. 

Numerical simulations have been performed to check the consistency of these results 
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Figure 12. CH simulations at fixed De. the drag coefficient / is found to increase with the 
Reynolds number. The comparison between the various curves shows that the drag coefficient 
reduces with the Deborah number. 



and their outcome is summarized in figure IT2l Here, the Deborah number has been fixed 
at different values and the drag coefficient has been plotted versus the Reynolds number. 
While / increases with Re, as expected, larger Deborah numbers are always found to be 
associated to smaller drag coefficients. 



9. Conclusions 

The weakly nonlinear dynamics of a viscoelastic Kolmogorov flow has been analysed 
both analytically and numerically. The physical reasons for considering this flow are that, 
despite the fundamental difference consisting in the absence of material boundaries, it 
has several analogies with channel flows and is one of the few well-known exact solutions 
of the Oldroyd-B model. 

The linear stability analysis for the Kolmogorov flow had already been developed by 
Boffett a et a!\ (|2005a) . No insights had however been given for the weakly nonlinear stage 
of evolution. This regime amounts to considering values of the Reynolds number close to 
the marginal stability curve separating stable from unstable regions of the phase-space. 
In the general nonlinear case (i.e. for arbitrarily large distances from the marginal curve), 
there is no way to solve the fully nonlinear equations. Conversely, close to the marginal 
curve, asymptotic perturbation techniques can be employed to capture the weakly non- 
linear dynamics. 

We found that the weakly nonlinear dynamics is described by Cahn-Hilliard-like equa- 
tions, with coefficients dependent on the Deborah number. The behaviour of these coeffi- 
cients with respect to De reveals that there exists a critical value of the Deborah number, 
where the system bifurcates to another regime. The resulting nonlinear equation still has 
a Cahn-Hilliard form, but contains a novel, fifth-order nonlinearity. 

Above the critical De, the "hydrodynamic" kink-antikink structures are replaced by 
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generalized structures. We have shown that their processes of annihilation are slowed 
down with respect to the standard Cahn-Hilliard equation. We also found a purely non- 
linear, subcritical mechanism of instability, which occurs for sufficiently large amplitudes 
of the initial perturbation. 

Our results demonstrate that, for hydrodynamical systems governed by a standard 
Cahn-Hilliard equation, the presence of an additional parameter might lead to higher- 
order nonlinear dynamics. A system where similar phenomena are to b e expected is the 
stratified Kolmogorov flow investigated bv lBalmforth &: Yound l|2005r) . with the role of 
elasticity played by stratification. 

Our results have been obtained both exploiting the multiscale expansion and via direct 
numerical simulations of the original equations and their coarse-grained version. The 
agreement between the Cahn-Hilliard dynamics and the full-resolved one is excellent 
even at large times. This is true for both the standard Cahn-Hilliard and the generalized 
one. The asymptotic analysis is thus able to capture all of the relevant features of the 
flow. 

In the last part of the work, we have presented some consequences for the problem 
of drag reduction. Although it is not common to talk about this effect in non-turbulent 
flows, we have shown that, even in the weakly nonlinear case, the injection of polymers 
induces a reduction in the drag coefficient, via the stabilization of the basic flow. Using 
the results of the nonlinear analysis, we have been able to give an analytical expression 
for the flow enhancement due to the polymers. The main qualitative conclusion is that 
drag reduction stems from a stabilization of the flow and appears to be a phenomenon 
coupling large and small scales. 
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Appendix A. Squire's theorem for Oldroyd-B 

Consider a parallel flow U = (U (z), 0). To investigate its stability properties one writes 
the linearized, nondimensional equations 

d t w + (u ■ V)w + (w ■ V)u = -Vg + /3Re~ 1 Aw + 

+ (1 - f3) Re' 1 De- 1 V • C (Al) 

8tC + (u • V)C + (w • V)<r = (Vuf ■ C + (Vw) T • a + 

+C ■ ( Vu) + <x • (Vw) - De^C (A 2) 

where w is the perturbation of the basic profile u, and £ is the perturbation of the basic 
stress tensor cr. 

We now perform a Fourier transform in the directions x and y, and in time, 

Wi(x,y,z,t) = j dLjdk x dk y e-^ t+kxX+kyy w t {k x ,k y ,bj, z) (A3) 

bj(x,y,z,t) = j ' du dk x dk v e- wt+k * x+k y y Qj (k x , k y , lj , z) (A4) 
Introducing the notation 



Cxz j - ( Cxx Cxy \ ^ ( &xz \ g ( ®xx O ' X y 

Cyz J ~ V Cyx Cyy J \ ^yz J \ °yx <Jyy 



the linearized equations in normal modes take the form 

du d? 
(-ioj + i k T ■ u)w + w z -^- = -ikq + f3Rer 1 (-k 2 + --^)w + 
dz dz z 



(A 5) 



+ (1 - /?) Re' 1 De' 1 (iz T • k + ^-tj (A 6) 



{-iu + ik T ■ u)w z = -j- + (3Re 1 {-k 2 + j^)w z + 

+ (1 - f3) Re- 1 De- 1 (V ■ t + ^C«) (A 7) 

, . ., t ^ -h. - d - du T du „ T 
(-jw + ik J • u + De L )z + w z — s = t • — + — • t 1 + 

dz dz dz 

+i(s ■ k)w T + iw(k T ■ s) + r-^-w T + ^-r T (A 8) 

dz dz 

(-iw + ik 2 -u + De t + w z — r = C zz ~r + 

ax dz 

+z(s • k)w z + iw(r T • k) + r-^-w z + -^-w (A 9) 

az az 

(-iw + ik T • u + De- X % z = 2i(r T ■ k)w z + 2^-w z (A 10) 

az 
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Consider the following transformation 



k x = \k\ w x = K^ w z = w z q = ^q 



Re=^Re De=^De ZJ = M w 



K k J -t 



|k| |k| ^ - | k | | k |2 S Z Z - | k |' 

From one can derive the equations for the overlined variables 

_ dU - 1 -2 d 2 

[-iuJ+ik x U(z)] w x + w z — = -i k x q + (3Re (-k x + "T - ^)^ ' 

+ (l- ^Tie^We^ (ik x ( xx + — \ 



(All) 



(A 12) 



[-ioJ + i k x U(z)] w z = -— + (3Re {-k x + -j^)w z + 



dz 2 ' 



-iu> + ik x ll(z) + De 



( I Rt 1 De 1 (ik x t x • -^C: 



— ds xx - dU — _ _ dw x 

C ra + w z — - — = 2t x — h 2ik x s xx w x + 2r x —-— 

dz dz dz 



(A 13) 
(A 14) 



-iuj + ik x U(z) + De 



- _ d _ - dU ._ - _ - 

t x + w z —r x — Q zz — h is xx k x w z + ik x w x r x 

dz dz 



dw z dw x 
+r x — 1 — - — (A 15) 

dz dz 



-iu> + ik x U(z) + De 



where we introduced the quantities 



- — n dw z 

Q zz = 2ik x r x w z + 2— — 
az 



k l-l± = l + -D-e 2 [U'{z)] 2 , r x 



^-L=De*U'(z) 



(A 16) 



(A 17) 



Equations IjA 12(1 - I|A 161) are exactly the same as (IA 6|) - (|A 10|) but with k y = 0, w y = 
0,Cxy = Cyy — Cyz = 0. Therefore they describe a two-dimensional linear disturbance 
of the basic flow at smaller Reynolds and Deborah numbers. If the three-dimensional 
perturbation w,£ is unstable at (Re,De), then the two-dimensional disturbance w,C is 
unstable at (Re, De) and its rate of growth is larger (Im(a7) > Im(w) > 0). 



Appendix B. Painleve analysis 

We perform a Painleve analysis to ascertain whether the fifth-order equation l|5.1|l is 
integrable as the usual cubic Cahn-Hilliard equation (|4.6I) . 

After rescaling dependent and independent variables, the stationary equation takes the 
form: 

-u-^--Xd 2 x u+lu 5 = 0. (Bl) 

The Painleve test consists in checking whether the structure of the solution around 
singularities in the complex plane has the form of a Laurent series. A simple balance of 
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the last two terms in the equation indicates that the singularity has order —1/2. The 
putative Laurent series should then be sought as: 



u(z) 



-1/2 



U0+U1Z + U2Z + u 3 z- + 



(B2) 



where z is the complex variable denoting the separation from the singularity z*. When 
the series (|B 211 is inserted into equation l|B a hierarchy of equations of the form 
a k u k — bk is obtained. The a^'s and b^s can be calculated in terms of Uk-i, ■ ■ ■ ,u . 
The impossibility for an arbitrary equation to have a Laurent series expansion is due to 
resonances, i.e. values of k such that = 0. Integrability is equivalent to checking that 
6fe = for the orders corresponding to resonances. In our case, it is easy to check that 

1/4 



-A k 



k 



u 



15A 

47 



a k = -A(fc+l)(fc-3) . (B3) 



The resonance is therefore at the third order and we need to perform the explicit calcu- 
lation up to that order to check whether or not 63 = 0. The algebra is elementary and 
the coefficients are: 



Ul = 



"2 = 



MO 

12A' A 
Using these values, one can verify that 

63 = 27^0^1 + 47UqUiU 2 
vanishes and the Painlcvc test is satisfied. 



1 

3 + 

Ul 



1287 



ulu 2 



Uqu\ 



(B4) 



(B5) 
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